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1.0  SUMMARY 


This  report  describes  the  application  of  large-eddy  simulation  /  Reynolds-averaged 
Navier-Stokes  (LES/RANS)  techniques  to  simulate  supersonic  flow  over  a  cavity  configuration 
tested  in  the  Air  Force  Research  Laboratory’s  (AFRL)  Trisonic  tunnel.  The  model  is  used  to 
compute  the  wavefront  aberrations  in  an  optical  beam  passing  through  the  cavity.  The  turbulence 
model  blends  a  RANS-type  closure  near  solid  walls  with  a  subgrid  model  in  the  free-stream 
based  on  the  ratios  of  estimated  inner  and  outer  turbulent  length  scales.  The  cavity  geometry  is 
modeled  using  an  immersed  boundary  method,  and  an  auxiliary  flat  plate  simulation  is 
performed  to  replicate  the  effects  of  the  wind-tunnel  boundary  layer  on  the  computed  optical 
path  difference.  Two-dimensional  proper  orthogonal  decomposition  modes  of  the  optical 
wavefront  are  computed  and  compare  favorably  with  wind  tunnel  data  despite  uncertainties 
about  inflow  turbulence  levels  and  boundary  layer  thicknesses  over  the  wind  tunnel  window. 
Dynamic  mode  decomposition  of  a  planar  wavefront  spanning  the  entire  cavity  reveals  that 
wavefront  distortions  are  driven  by  shear  layer  oscillation  at  the  Rossiter  frequencies;  these 
disturbances  create  eddy  shocklets  that  propagate  into  the  free-stream  and  create  additional 
optical  path  disturbances. 

It  should  be  noted  that  the  subject  discussed  in  this  report  is  different  from  that  originally 
proposed,  which  focused  on  simulations  of  the  effects  of  heat  release  due  to  lasing  chemistry  on 
shock-train  formation  within  chemical  oxygen  iodine  lasers  (COILs).  A  delay  in  funding  (due  to 
the  DoD  sequester)  led  to  a  change  in  the  scope  of  work,  spurred  by  Ilya  Zilberter’s  stay  at 
AFRL  during  the  summers  of  2014  and  2015.  Ilya  was  the  Ph.D.  student  supported  under  this 
work.  A  no-cost  extension  has  been  awarded  for  this  work  to  allow  us  to  complete  some  of  the 
work  proposed  under  the  original  statement  of  work. 

2.0  INTRODUCTION 

Externally-mounted  optical  systems  (e.g.  an  aircraft-mounted  laser)  may  experience  disturbances 
in  the  beam  path  caused  by  turbulence  or  other  flow  features.  In  the  case  of  an  optical  device 
inside  an  open  cavity,  the  beam  must  pass  through  a  turbulent  detached  shear  layer  and  a  series 


Aperture  plane 

Figure  1.  Schematic  of  cavity  flow  with  optical 
reflector  at  bottom  of  cavity 

of  acoustic  waves  which  result  from  flow  re-attachment  at  the  back  of  the  cavity,  both  of  which 
can  affect  the  coherence  and  power  of  the  beam  by  aberrating  the  beam  wavefront  (Figure  1) 
When  designing  such  systems,  it  is  important  to  be  able  to  predict  and  model  these  aero-optical 
effects. 
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The  present  study  aims  to  do  so  using  a  hybrid  Large-Eddy  Simulation  /  Reynolds- 
Averaged  Navier-Stokes  (LES/RANS)  turbulence  closure.  This  approach  directly  resolves  large 
turbulent  eddies  in  the  free  stream  while  switching  over  to  a  RANS  model  to  simulate  the 
smaller-scale  eddies  near  the  walls.  Thus  the  simulation  should  be  able  to  capture  the  transient, 
large-scale  turbulence  structures  characteristic  of  these  types  of  cavity  flows  while  requiring 
smaller  computational  meshes  than  traditional  LES.  Detached  eddy  simulation  (DES)  and 
LES/RANS  methods  have  been  applied  to  cavity  flows  in  several  prior  studies,  including  for 
both  supersonic  and  subsonic  cavity  flows,  and  are  generally  able  to  capture  the  cavity  flow 
physics  with  good  accuracy.  Among  more  recent  efforts,  Wang,  et  al.  [1]  used  a  hybrid 
LES/RANS  model  to  simulate  a  relatively  high  aspect-ratio  cavity  in  Mach  1.8  and  2.5  free- 
stream  conditions  and  were  able  to  accurately  predict  the  oscillation  frequencies  of  the  shear 
layer  as  well  as  mode  switching  to  more  wake-like  behavior.  Peng  and  Liecher  [2]  used  a  hybrid 
LES/RANS  code  to  simulate  a  three-dimensional  Mach  0.85  cavity  flow  including  the  entire 
wind  tunnel  geometry.  While  their  approach  managed  to  capture  the  first  three  acoustic 
disturbance  modes,  the  spectral  behavior  associated  with  higher  frequencies  was  not  sufficiently 
resolved,  partly  due  to  insufficient  resolution  of  the  computational  grid  and  simulation  time  step. 

An  advantage  of  LES  and  LES/RANS  methods  is  that  the  direct  simulation  of  large  eddy 
structures  allows  for  the  prediction  of  aero-optical  effects,  which  arise  from  density  fluctuations 
within  the  optical  path.  Mani,  et  al.  [3]  argue  that  the  smallest  optically-active  turbulence  scales 
fall  within  the  range  typically  resolved  in  standard  LES  applications.  [4]  Morgan  and  Visbal  [5] 
used  a  RANS/  (implicit)  LES  turbulence  model  to  simulate  aero-optical  aberrations  resulting 
from  flow  over  an  optical  turret,  achieving  fair  comparison  with  experimental  data  over  a  range 
of  elevation  angles,  albeit  with  deviations  in  the  separated  flow  region  behind  the  turret.  Engfer, 
et  al.  [6]  performed  delayed  detached-eddy  simulations  (DDES)  of  flow  around  the  aircraft- 
mounted  SOFIA  telescope  but  their  results  were  hampered  by  insufficient  turbulence  resolution 
in  the  inflow  boundary  layer,  which  caused  computed  wavefront  aberrations  to  differ  from 
experimental  values  by  a  factor  of  4.  Additionally,  the  instrument  used  to  detect  spatially- 
resolved  wavefront  aberrations  in  the  experiment  was  out  of  focus  and  failed  to  produce  evidence 
of  spatially  coherent  turbulent  structures. 

Presently,  we  examine  a  supersonic,  three-dimensional  cavity  flow  on  a  high-resolution 
computational  mesh  in  order  to  test  the  capability  of  LES/RANS  models  to  predict  aero-optical 
phenomena.  While  most  prior  studies  have  used  experimental  validation  in  the  form  of  pressure 
data  sampled  along  the  cavity  floor  and  walls  or  one-dimensional  wavefront  deflection  spectra, 
recent  experiments  in  the  Trisonic  Gasdynamic  Facility  (TGF)  [7]  at  Wright-Patterson  AFB  have 
yielded  Shack-Hartmann  wavefront  measurements  for  a  range  of  supersonic  cavity  flows  that 
directly  show  the  optical  aberrations  caused  by  the  detached  shear  layer.  Using  these  data  for 
validation  allows  for  the  comparison  of  not  only  the  modal  frequencies  of  the  optical 
disturbances,  but  of  their  shape  as  well.  Figure  2  shows  the  cavity  geometry  within  the  blocked- 
out  simulation  domain.  The  optical  beam  is  highlighted  in  red. 
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Figure  2.  Overview  of  simulation  domain  with  optical  beam  path  highlighted  in  red 


3.0  METHODS,  ASSUMPTIONS,  AND  PROCEDURES 
3.1  Time  Advancement 

NCSU’s  REACTMB  flow  solver  is  used  in  the  present  effort.  REACTMB  solves  the 
Navier-Stokes  equations  governing  a  multi-component  mixture  of  gases  on  simply-connected, 
multi-block  structured  meshes  using  finite-volume  methods.  REACTMB  is  designed  for  use  on 
massively-parallel  computers  and  uses  MPI  for  message-passing.  Large  structured  meshes  are 
decomposed  into  a  number  of  smaller  blocks,  which  are  then  partitioned  over  the  number  of 
requested  cores  using  a  simple  divide-and-conquer  strategy.  REACTMB  discretizes  the  Navier- 
Stokes  system  in  time  using  a  Crank-Nicholson  approach: 

Tjn+\,k  _Tjn  i  i 

Q- - — +  -(l  +  0)R(Vn+lk )  +  -(!- 0)R(Vn)  =  O  (1) 

At  2  2 

where  Q  is  the  cell  volume,  At  is  the  time  step.  U  is  the  vector  of  conserved  variables,  and  R  is 
the  residual  vector.  The  function  0  is  defined  as 

6>  =  -(l-tanh(^^)),  dt  =lxl0‘4m,  Ad  =  0.2dt  (2) 

2  Ad 

Here,  d  is  the  distance  to  the  nearest  solid  surface.  The  function  6  switches  the  time 
discretization  from  Crank-Nicholson  to  Euler  implicit  for  mesh  cells  essentially  within  the 
laminar  sub-layer.  Some  loss  of  temporal  accuracy  results,  but  this  approach  is  necessary  to 
suppress  oscillations  in  the  pressure  and  transverse-velocity  fields  for  mesh  cells  with  a  very  high 
aspect  ratio.  The  matrix  system  that  results  from  linearizing  (1)  is  solved  approximately  using 
a  block  incomplete  lower-upper  (ILU)  factorization  method,  and  the  system  is  converged  to  a 
prescribed  tolerance  over  a  sequence  of  sub-iterations.  Jacobian  matrix  elements  are  stored  over 
the  number  of  blocks  mapped  to  a  particular  processor,  allowing  the  “freezing”  of  the  matrix 
elements  and  their  factorization  over  the  duration  of  the  sub-iterations.  This  reduces  the 
computational  workload  significantly. 
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3.2  Flux  Formulation 


To  sustain  turbulence,  it  is  necessary  to  reduce  numerical  dissipation  significantly.  The 
strategy  employed  in  REACTMB  combines  a  TVD  scheme  with  a  fourth-order  central  difference 
scheme.  Edwards’s  low  diffusion  flux-splitting  scheme  (LDFSS)  [8]  is  used  as  the  Riemann 

solver.  The  primitive-variable  vector  W  =  [ps,u,v, w,T,k,o)]T is  used  in  the  reconstruction. 

TVD  schemes  preserve  monotonicity,  which  compromises  their  ability  to  resolve  small-scale 
turbulence.  One  means  of  alleviating  this  problem  is  to  blend  the  monotonicity-preserving  left- 
and  right-state  values  with  an  averaging  operator  that  leads  to  a  fourth-order  central  difference, 
so  that  the  latter  is  used  in  regions  of  high  vorticity  (boundary  layers,  shear  layers)  and  the 
former  is  used  in  more  ‘inviscid’  regions,  where  strong  shocks  might  be  present.  A  function 
due  to  Ducros,  et  al.  [9],  defined  at  a  mesh  cell  as 

f= - . - T,  f^lxlO-8^ /max(Ax,Ay,  Az)  ,  (3) 

(V-V)2  +  \(o\2  +s2 


is  used  for  this  purpose.  Here,  the  divergence  of  velocity  is  compared  with  the  vorticity  value.  If 
the  latter  is  much  larger  (in  shear  and  boundary  layers,  for  example),  the  function  moves  toward 
zero,  and  in  free-stream  regions  near  shocks,  the  function  approaches  one.  At  a  particular  cell 
interface  i+1/2,  we  use  the  function  as  follows: 


”Vi+ 1/2  —  "R, i+1/2  + 

wA  =WA 

yyL,i+ 1/2  yy  R,  i+1/2 


■- — (w  +w;.+1)-— ( 

12  '  '  12 


(4) 


The  superscripts  M  and  A  refer  to  the  monotonicity-preserving  scheme  and  the  fourth-order 
averaging  operator,  respectively.  This  scheme,  denoted  as  LD-TVD  for  low-dissipation  TVD,  is 
used  for  all  LES/RANS  calculations  presented  in  this  paper.  Viscous  and  diffusive  terms 
appearing  in  the  equation  system  are  discretized  using  second-order  central  differences. 


3.3  LES/RANS  Model 


NCSU’s  hybrid  LES/RANS  methodology  [10]  is  used  as  the  baseline  for  the  current 
work.  In  the  LES/RANS  model,  the  effects  of  anisotropic  near- wall  eddies  are  modeled  using 
RANS  concepts  (Menter’sk-&> baseline  model),  whereas  the  larger  turbulent  eddies  away  from 
solid  surfaces  are  captured  using  a  large-eddy  simulation  method.  The  shift  between  the  closure 
models  is  facilitated  by  modifying  the  eddy  viscosity  field  according  to 


k 

1 

r  i  ^ 

A,  =  P 

a-ix^+r- 

CO 

•r-2 

t-H 

1 

,  |«5! 

VI 

1 

Here,  T  is  a  time-dependent  blending  function  that  connects  the  RANS  and  LES  branches  and  the 
quantity  XN  is  a  ratio  of  outer-  to  inner-layer  turbulence  length  scales,  with  the  former  calculated 

using  both  ensemble-averaged  and  instantaneous  turbulence  data.  Specific  definitions  are  as 
follows: 

K  =g(L,JLjL„„  *(/„„„)  =  min[D„max(l,D2^  -)]  (6) 

l  ,  V  co 

outer  1 

where 
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\\Ovco  +  k+k^ 

‘'outer  '-'N,  ^-,1/2 —  ’  inner 

y  coco 

The  model  constants  D{  and  D2  are  assigned  values  of  10  and  0.5,  based  on  calibrations  for  flat- 
plate  boundary  layers.  The  mesh  scale  is  taken  to  be  the  maximum  spacing  over  all  three 

coordinate  directions.  This  form  serves  to  shift  the  closure  to  unsteady  RANS  when  there  is  no 
possibility  of  resolving  the  largest  turbulence  length  scales.  If  the  maximum  mesh  scale  is 
selected  as  the  outer-layer  scale  consistently,  then  the  LES/RANS  model  behaves  similarly  to 
detached-eddy  simulation,  serving  primarily  to  isolate  turbulent  boundary  layers  from  massively- 
separated  regions.  The  quantity  kR  is  the  resolved  turbulence  kinetic  energy,  calculated  via 
ensemble-averaging  as 

PK  =  \  (pukuk  ~  pUkfUk )  (7) 

2  p 

The  quantities  k  and  co  are  ensemble-averages  of  the  modeled  turbulence  kinetic  energy  and 
specific  dissipation  rate  variables,  which  are  obtained  from  Menter’s  model.  The  use  of 
ensemble-averaged  as  well  as  instantaneous  data  allows  the  RANS-to-LES  transition  location  to 
‘float’  about  a  time-mean  location  that  is  a  function  of  the  ensemble-averaged  state  of  the 
boundary  layer.  This  enables  the  blending  function  to  respond  more  directly  to  large  changes  in 
the  turbulence  length  scales  that  can  result  from  shock  interactions.  Ensemble-averages  are 
currently  calculated  using  an  exponentially- weighted  moving  average:  Qn  =  (l- A)Q "_1  +  AQ" 
with  A  =  At/r.  The  time  scale  r  is  defined  as  follows  for  the  cases  considered  herein: 


T  =  nan(t,tres),  t<4t, 
=  f-3r,r>4/ 


(8) 


where  rres  =  LI u.r  is  defined  in  terms  of  the  length  of  the  domain  L  and  the  free-stream  velocity 
ux  This  form  assumes  that  a  statistically-stationary  state  will  emerge  after  about  four  residence 
times 


3.4  Immersed  Boundary  Method 

Accurate  and  affordable  LES  simulations  require  fine,  isotropic  mesh  resolution  in 
regions  of  interest.  It  is  desirable  to  avoid  mesh  clustering  that  extends  away  from  solid  surfaces 
into  ‘free-stream’  or  ‘ffee-shear  layer’  areas  -  the  cells  that  result  are  normally  anisotropic  and 
can  lead  to  severe  time-step  restrictions  as  well  as  oscillations  where  the  fourth-order  central 
different  method  is  used.  The  ‘boat-shaped’  cavity  geometry  shown  in  Figure  2  has  several 
sharp  edges  where  different  flat  surfaces  are  joined  -  resolving  these  with  typical  structured 
meshes  leads  to  the  problems  mentioned  above,  which  greatly  degrade  the  accuracy  and 
efficiency  of  the  flow  solver.  Attempts  were  made  to  generate  body-fitted  meshes  using  both 
GridPro  and  GridGen,  but  a  successful  rendering  could  not  be  obtained.  Instead,  we  opted  to 
utilize  NCSU’s  established  immersed-boundary  method  [11-13]  to  mimic  the  effects  of  the 
object  on  the  surrounding  flow.  A  220  M  cell  Cartesian  mesh  containing  isotropic  cells  of  sizes 
no  greater  than  0.2  mm  in  the  region  of  interest  (the  boat-tail  geometry  and  its  cavity)  was 
generated,  the  boat-tail  geometry  rendered  as  a  stereolithography  (STL)  file,  and  procedures 
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described  in  [11]  and  [12]  used  to  ‘embed’  the  object  into  the  computational  domain.  The 
‘embedding’  process  involves  first  the  calculation  of  the  signed  distance  from  any  mesh  cell  to 
the  nearest  point  on  the  STL  rendering.  Cells  inside  the  object  are  classified  as  ‘interior’  cells  - 
fixed  properties  are  enforced  for  these  cells  and  they  do  not  interact  with  other  cells  in  the 
domain.  Cells  outside  the  object  but  with  no  face,  edge,  or  vertex  neighbor  inside  the  object  are 
classified  as  ‘field’  cells  -  the  compressible  Navier-Stokes  equations  are  solved  in  these  cells. 
Cells  outside  the  object  but  with  at  least  one  face,  edge,  or  vertex  neighbor  inside  the  object  are 
called  ‘band  cells’  -  here,  flow  properties  are  reconstructed  based  on  the  external  field-cell 
solution  and  on  the  prescribed  boundary  conditions  at  the  embedded  surface. 

The  following  first-order  accurate  closures  are  used  for  the  fluid  properties  in  the  band 
cells,  where  the  subscript  ‘I’  indicates  properties  obtained  at  an  interpolation  point  located  along 
the  normal  line  extending  outward  from  the  nearest  surface  location  corresponding  to  the  band 
cell  in  question,  and  the  subscript  “B”  indicates  the  band-cell. 

Pb  =  p{dj) 

I'd  \k 

uB,i~us,i=uT,i(di)  ~r  +uNi(dI)c(p,dj,dB), 

d,  J  (9) 

unMi^>  =  (uj(dl)~uS,j)njni’ 
uT,i  (dI )  =  (ui  (dI )  -  US ,i )  -  u N,i  (dI  > 

In  these  expressions,  n,  is  the  normal  vector  at  the  closest  point  on  the  body  surface,  d  is  a 
distance  from  the  nearest  surface  point,  us  •  is  the  velocity  at  the  nearest  surface  point,  and  k  is  a 

power-law.  The  choice  of  k  allows  the  model  to  replicate  a  turbulent  velocity  profile  (k=  l/7  or 
1/9)  or  a  laminar  profile  (k=  1).  To  obtain  the  temperature  distribution  near  the  surface,  Walz’s 
relation  for  the  temperature  distribution  within  a  compressible  boundary  layer  is  used: 

Isothermal  wall 


Adiabatic  wall: 


In  this,  r  is  the  recovery  factor  and  [uT  i{dj)f  is  the  kinetic  energy  associated  with  the  tangential 

velocity  component  at  the  interpolation  point.  The  function  c{p,  dl  ,dB)  that  scales  the  normal 

velocity  component  in  Eq.  (9)  is  determined  by  enforcing  a  discrete  form  of  the  continuity 
equation  at  each  band  cell  using  a  locally-parallel  flow  assumption.  Details  are  given  in  Ghosh 
etal.  [13] 

The  turbulence  variables  in  the  band  cells  are  defined  as 


Tm  _  ,  L  ,  r(r- 1)  (d)f\\^L 

T(d, )  T{dj)  {  T(dj)  2  yRT(dj)  T’1  1  ){dj 


[Mr  i  idi  )f\  — 
2rRT(d{)  T’1  (d, 


-1b—  =  \+  rSl  [U  (dj)]2 (1 -[  —  }  ) 
T(dj)  2  yRTidj)  \dI  J 
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kB=u^/^jc^  ,  ojb=  uT  K^C^k<1b)  :  d+  >10.934 


kB  =  k(dj ) 


\2 


ydj  J 


,  coB  =  60vw/(0.015dl) :  d+  <  10.934 


d+  =uTdb/vw,  uT  =\uTi(dI)\/(ln(d+)/ic  +  5A) 
(iterative  solution) 


(12) 


To  arrive  at  this  form,  we  assume  equivalence  between  the  result  provided  by  the  power-law 
profile  and  the  law  of  the  wall  within  the  band  cells.  Details  regarding  the  calculation  of  flow 
properties  at  the  ‘interpolation  point’  may  be  found  in  [11-13]. 


3.5  Optical  Path  Difference  Calculations 

Optical  Path  Difference  (OPD)  data  from  the  physical  experiment  at  the  TGF  is  provided 
by  a  3cm  diameter  laser  beam  which  is  reflected  off  of  the  bottom  of  the  cavity  and  into  a  Shack- 
Hartmann  wavefront  detector  outside  of  the  wind  tunnel.  The  OPD  at  each  point  on  the  detector 
is  related  to  the  Optical  Path  Length  (OPL), 


OPEHA,  z,  t )  =  OPUx,  z,  t )  -  OPldX  z,t)  (13) 

which  in  turn  is  proportional  to  the  integral  of  the  gas  density  in  the  beam  path, 

b 

OPL(x,  z,  t)  \kgdp(x,  y,  z,  t)  +  l]dz  (14) 


where  p  is  the  density  and  Kgd  is  the  Gladstone-Dale  constant  (2.27xl0'4  m3/kg  for  visible  light 
in  air).  Time-accurate  density  data  from  the  LES/RANS  simulations  can  be  used  to  reconstruct 
the  simulated  OPD  along  a  virtual  beam  path  and  directly  compared  with  sensor  data  from  the 
physical  experiment.  To  this  end,  path  integrals  of  the  density  along  a  3 -cm  wide  region  above 
the  center  of  the  simulated  cavity  are  sampled  at  a  rate  of  24  kHz  as  in  the  physical  experiment. 
Each  integral  is  then  processed  to  convert  the  density  into  an  OPD  snapshot  according  to  Eqs.  13 
and  14.  The  OPD  is  further  processed  to  remove  the  tip  and  tilt  components  of  the  wavefront 
which  represent  the  average  gradients  of  OPD  in  the  two  directions  across  the  aperture. 
Traditionally,  this  is  done  in  experimental  studies  to  remove  the  effects  caused  by  vibrations  in 
the  wind  tunnel  or  optical  instruments,  which  would  cause  mono-directional  deflections  in  the 
measured  wavefront.  However,  it  is  understood  that  some  flow  features  introduce  tip  and  tilt 
components  that  do  not  fall  under  experimental  error  [7].  For  example,  the  thickening  of  a 
boundary  layer  across  the  aperture  causes  a  significant  tilt  in  the  wavefront.  Nonetheless,  even 
though  OPD  data  from  the  LES/RANS  simulations  is  free  of  instrumentation  error,  the  tip  and 
tilt  components  are  still  filtered  out  for  the  sake  of  comparison  with  experimental  results.  This  is 
done  by  minimizing  a  target  function  G  over  the  aperture  Ap  [4]: 

G  =  \\\pVlAx,z,t)-(Ax  +  Bz  +  C)2\lxdz  (15) 

Ap 
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x(m) 

Figure  3.  Superimposed  LES/RANS  centerline  snapshots  of  cavity  and  wind  tunnel  wall  domains. 

Optical  path  outlined  in  white 

Solving  dG/d(A,B,C)  =  Ofits  a  linear  OPL  wavefront  to  the  actual  data  in  a  least-squares 
sense.  The  tip/tilt-removed  OPD  is  then  given  by: 

OPD(X  z,t)  =  OPD(x,  z,  t)  -  ( Ax  +  Bz  +  C)  (16) 

where  (Ax+Bz  +  C)  represents  a  linear  OPL  function  and  implicitly  includes  the  aperture- 
averaged  OPL  term  in  Eq.  5.  In  a  sense,  this  tip/tilt  removal  procedure  is  a  basic  form  of  adaptive 
optics  that  results  in  a  zero-mean  OPD  wavefront  that  only  includes  higher-order  aberrations. 

Because  the  top  wall  of  the  wind  tunnel  is  not  modeled  in  the  cavity  simulation,  a 
separate  turbulent  boundary  layer  simulation  with  identical  free-stream  conditions  was 
performed  on  a  4  million  cell  mesh  to  account  for  the  effect  of  turbulence  along  the  wind  tunnel 
window  on  the  beam.  A  superimposition  of  simulation  snapshots  along  both  the  cavity  and 
boundary  layer  centerlines  is  shown  in  Figure  3,  with  the  Shack-Hartmann  beam  location 
outlined  in  white  to  highlight  the  turbulent  regions  that  the  beam  passes  through.  The  boundary 
layer  thickness  of  18mm  and  rms.  OPD  levels  are  compared  to  empty  wind  tunnel  data  given  in 
Ref.  [7]  for  validation,  and  the  density  integral  along  the  beam  path  is  added  to  that  of  the  cavity 
simulation  when  computing  the  total  OPD. 

4.0  RESULTS  AND  DISCUSSION 

4.1  Flow  Structure 

The  LES/RANS  cavity  and  boundary  layer  flow  simulation  were  initialized  from  a 
converged  RANS  solution  and  progressed  for  approximately  35  cavity  flow-over  times,  yielding 
both  instantaneous  and  time-averaged  flow  data  which  are  useful  in  giving  an  overview  of  the 
overall  flow  dynamics.  Figure  4  compares  a  Schlieren  photography  snapshot  from  the  physical 
experiment  with  instantaneous  density  gradient  (dp/dy)  contours  along  the  centerline  of  the 
simulation  domain.  Overall,  the  LES/RANS  captures  the  broad  structure  of  flow  over  the  body. 
The  angles  of  the  oblique  shocks  at  the  splitter  tip  and  cavity  edge  are  well-replicated,  as  is  the 
region  of  shed  shocklets  above  the  shear  layer  spanning  the  cavity.  These  are  caused  by  the  high¬ 
speed  free-stream  air  flow  impacting  slower-moving  eddies  as  they  rise  above  the  shear  layer. 
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Figure  4.  Schiieren  snapshot  of  domain.  Left:  experiment.  Right:  LES/RANS  centerline 


The  'wave'-like  shock  emanating  from  the  rear  cavity  wall  is  also  present  and  is  a  product  of  the 
unsteady  impingement  of  the  shear  layer  on  the  rear  wall.  As  desired,  the  geometry  of  the  body 
creates  a  turbulent  shear  layer  upstream  of  the  cavity  without  the  need  for  artificial  fluctuations. 
In  the  LES/RANS  snapshot,  the  oblique  shocks  above  the  body  appear  to  lose  strength  towards 
the  top  of  the  domain  because  of  numerical  dissipation  caused  by  larger  cell  sizes  towards  the  top 
boundary.  The  contraction  region  underneath  the  body  also  exhibits  a  noticeable  'stair-stepping' 
effect  due  to  the  IB  method  attempting  to  resolve  a  curved  surface  using  rectilinear  cells.  This 
effect  is  a  drawback  of  enforcing  full  mass  conservation  within  the  'band  cells'.  Although  this  is 
not  ideal,  the  flow  underneath  the  body  does  not  affect  conditions  near  the  cavity  itself  and  so  the 
effect  was  not  considered  important. 

Figure  5  shows  the  instantaneous  density  and  stream-wise  velocity  contours  at  the  cavity 
centerline.  The  shear  layer  rapidly  thickens  and  develops  very  fine  turbulent  structure  but  has 
enough  energy  to  span  the  cavity  without  significant  deflection.  Based  on  the  high  level  of  small- 
scale  density  fluctuations  in  the  shear  layer,  it  can  already  be  expected  that  this  region  of  the  flow 


Figure  5.  Centerline  snapshot  of  cavity  from  LES/RANS  simulation.  Left:  density  contours.  Right: 

streamwise  velocity  contours 


will  produce  the  most  aero-optical  aberrations  in  the  path  of  the  beam.  Significant,  but  larger- 
scale,  density  variations  are  also  seen  above  the  shear  layer  in  the  region  of  shed  shocks.  In  the 
density  snapshot,  a  density  wave  can  be  observed  inside  the  cavity  at  x=0.25m.  As  the  wave 
travels  upstream  and  hits  the  front  cavity  wall,  it  will  create  a  deflection  in  the  shear  layer  that 
propagates  the  shear  layer  oscillation  cycle. 

The  centerline  velocity  contour  (Figure  5)  reveals  a  strong  region  of  reversed  flow  in  the 
downstream  half  of  the  cavity,  especially  at  the  aft  wall  and  along  the  cavity  floor.  The  flow  in 
the  cavity  is  actually  divided  into  four  distinct  recirculation  zones  which  are  symmetric  about  the 
centerline.  The  strong  reversed  flow  at  the  aft  wall  caused  by  the  shear  layer  impingement  moves 
along  the  cavity  bottom  and  feeds  into  the  two  vortices  in  the  upstream  corners;  the  air  in  these 
zones  then  rises  until  it  is  entrained  by  the  shear  layer  and  carried  either  out  of  the  cavity  or  into 
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the  two  downstream  vortices.  The  complexity  of  this  flow  structure  compared  to  the  single 
recirculation  zones  normally  predicted  in  two-dimensional  cavity  simulations  highlights  the  need 
to  consider  this  class  of  problems  as  fully  three-dimensional. 

Finally,  Figure  6  shows  an  instantaneous  density  isosurface  near  the  top  of  the  shear  layer 
to  help  visualize  the  some  of  the  turbulence  structure  in  the  optical  path  (shown  in  red).  As  the 
shear  layer  detaches  from  the  front  cavity  wall,  it  is  largely  two-dimensional  in  nature.  However, 
by  the  time  the  flow  encounters  the  beam  path  it  has  devolved  into  large,  uncorrelated,  three- 
dimensional  eddies.  The  largest  of  these  spans  about  half  of  the  beam  width;  the  aero-optical  data 
from  the  virtual  Shack-Hartmann  sensor  should  therefore  reflect  all  or  most  of  the  turbulent 
scales  in  the  flow,  although  not  necessarily  all  of  the  larger-scale  flow  dynamics. 


0  25 


0.35 


Figure  6:  Instantaneous  density  isosurface  above  cavity  colored  by  streamwise  velocity.  Optical 

path  shown  in  red. 


4.2  Cavity  Spectral  Behavior 

In  order  to  quantify  the  oscillations  of  the  shear  layer,  wall  pressures  along  the  cavity  line 
were  sampled  at  500kHz  and  processed  into  pressure  spectra.  Figure  7  shows  pressure  spectra  at 
the  center  of  the  front  and  rear  cavity  walls.  The  dotted  lines  denote  the  dominant  frequencies 
predicted  by  the  empirical  Rossiter  formula  [14] 
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Front  wall  pressure  power  spectral  density 


5000  6000 
Frequency  (Hz) 


Frequency  (Hz) 


Figure  7.  Power  spectral  density  at  center  of  front  and  rear  cavity  lip.  Dotted  lines  show  first  five 
oscillation  frequencies  predicted  by  Rossiter’s  formula. 


f  = ^2L 

Jm  L 


m  ■ 


■0.25 


M. 


>(l  +  + 


Vl/2 


1 

0.58 


(17) 


where  m  is  the  mode  number,  Uv  and  My  are  the  free-stream  velocity  and  Mach  number,  and  L  is 
the  cavity  length.  Table  1  lists  the  first  five  Rossiter  Table  h  First  flve  Rossiter  frequencies 


modes  for  this  cavity.  For  both  locations,  the  dominant  based  on  Eq.  17 

frequencies  closely  align  with  the  Rossiter  modes,  with  _ 

the  front  wall  displaying  first  mode  dominance  and  the  Mode  Frequency 

rear  wall  aligning  with  the  second  mode.  At  both  number _ (Hz) 

stations,  a  significant  low-frequency  component  is  1  816 

observed.  This  is  likely  due  to  transitions  to  'wake  mode'  2  1904 

behavior,  which  is  caused  by  the  boundary  layer  3  2993 

upstream  of  the  cavity  intermittently  detaching  and  4  4081 

shedding  a  large  vortex.  Wake  mode  behavior  has  been  _ 5 _ 5170 

observed  in  both  numerical  [15]  and  experimental  [16] 


studies  of  similar  cavities  at  supersonic  flow  conditions.  Overall,  the  spectral  data  agree  very 
well  with  both  the  empirical  oscillation  frequencies  and  the  behavior  seen  in  prior  studies. 

4.3  Aero-Optical  Analysis 

Density  data  from  cells  along  the  beam  was  sampled  every  208  iterations  (equivalent  to 
24  kHz,  which  is  the  sampling  rate  of  the  Shack-Hartmann  detector)  and  used  to  compute  the 
OPD  according  to  Eqns.  13-16.  In  order  to  include  the  effects  of  the  wind  tunnel  wall  in  the 
cavity  simulation,  the  density  integrals  from  the  boundary  layer  and  cavity  simulations  were 
linearly  added  before  computing  the  OPD.  The  resulting  path  length  was  greater  than  the 
distance  between  the  cavity  and  tunnel  wall  in  the  experiment;  however,  this  included  a  large 
laminar  free-stream  region  which  contributed  nothing  to  the  OPD  when  the  average  OPL  was 
subtracted  from  the  solution  as  in  equation  4.  In  effect,  the  total  length  of  the  beam  is  arbitrary 
when  computing  the  OPD  as  long  as  the  optically  active  regions  of  flow  (i.e,  regions  containing 
turbulent  eddies  or  shock  waves)  are  included  in  the  path.  The  cross-section  of  the  3cm  beam 
was  resolved  in  30x30cells  in  the  boundary  layer  simulation  and  120x120  cells  in  the  cavity 
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simulations.  When  computing  the  joint  OPD,  the  cavity  data  was  interpolated  onto  the  coarser 
grid.  This  was  still  sufficient  resolution  for  experimental  validation  as  the  Shack-Hartmann 
detector  resolved  the  beam  using  32x32  apertures.  Subsequent  analysis  of  OPD  from  the  cavity 
simulation  that  excludes  the  wall  boundary  layer  is  performed  on  snapshots  with  the  original 
120x120  resolution. 

A  comparison  of  the  OPD  data  from  the  LES/RANS  simulation  and  physical  experiment 
is  given  in  Figure  8.  Subfigure  (a)  shows  an  instantaneous  snapshot  of  the  tip  and  tilt  removed 
OPD  from  the  simulation.  The  time-averaged  mean  OPD  has  also  been  removed,  resulting  in  a 
snapshot  capturing  the  unsteady,  higher-order  aberrations.  Higher  values  of  OPD  indicate 
regions  of  the  beam  that  passed  through  higher  density  areas  and  experienced  a  relative  phase 
delay.  High  levels  of  disturbances  in  the  wavefront  can  be  observed  and  are  a  direct  reflection  of 
turbulent  eddies  in  the  beam  path.  Figures  8(b)  and  8(c)  compare  the  RMS  values  of  the 
mean/tip/tilt  removed  OPD,  calculated  as: 

OPDBMS.i|VOPDf-OPD_,  0PD_4|0PD,  US, 


where  N  is  the  total  number  of  OPD  snapshots.  Overall,  the  structure  of  OPD  variance  in  the 
aperture  is  well-replicated  in  the  simulation,  barring  a  region  of  higher  phase  delay  in  the  center 
of  the  aperture  in  the  experiment.  There  is  no  ready  explanation  for  the  discrepancy,  although  it 
could  be  due  to  wind  tunnel  effects.  The  computed  OPDrms  levels  are  also  as  much  as  20% 
locally  higher  than  the  experimental  values.  Commonly,  an  ensemble-averaged  value  of  OPDrms 
over  the  aperture  is  computed  to  provide  a  scalar  measure  of  optical  distortion: 


™  i= l 


i  M 

—  XtOPD^-OPD^)2 


1/2 


(19) 


where  M  is  the  number  of  sub-apertures  in  each  snapshot.  The  values  of  OPDrms, s  for  the  full 
optical  path,  as  well  as  the  isolated  wall  boundary  layer  and  cavity  simulations,  are  given  in 
Table  2.  The  levels  of  distortion  caused  by  the  tunnel  boundary  layer  are  half  that  of  the 
combined  flow,  although  still  significant.  The  computed  value  of  7.30  nm  for  the  boundary  layer 


x(m) 

(a)  (b)  (c) 

Figure  8.  OPD  levels  along  beam  path,  (a)  RMS  OPD  from  experiment,  (b)  RMS  OPD  from 
LES/RANS  simulation,  (c)  Snapshot  of  mean,  tip,  and  tilt  -removed  OPD  from  simulation 

is  slightly  lower  than  the  experimental  value  of  8.0  nm  presented  in  [7].  At  the  time  of  writing 
the  experimental  OPDrms, s  values  for  the  cavity  flow  were  not  available;  however,  due  to  the 
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overall  higher-than-expected  levels  of  OPDrms  seen  in  the  LES/RANS  data  in  Figure  8,  it  is 
expected  that  the  values  from  the  simulation  are  moderately  higher  than  those  of  the  experiment. 

Table  2.  Computed  OPDrms, s  values  over  aperture 


Domain 

OPDrms, s  (nm) 

Combined  cavity  and  wind  tunnel  wall 

14.19 

Wind  tunnel  wall  only 

7.30 

Cavity  only 

12.07 

4.4  Proper  Orthogonal  Decomposition  analysis  of  OPD 

To  better  understand  the  turbulent  structure  driving  the  optical  aberrations,  the  snapshots 
of  OPD  sampled  over  the  course  of  the  simulation  can  be  decomposed  into  proper  orthogonal 
(POD)  modes.  POD  analysis  has  widely  been  used  in  studies  of  both  aero-optics  [7,17]  as  well  as 
generalized  fluid  flows  [18].  Each  spatial  mode  is  a  basis  function  that  approximates  the 
variation  in  the  wavefronts  in  a  least-squares  sense.  Once  the  POD  for  a  flow  has  been 
calculated,  each  snapshot  can  be  reconstructed  using  a  linear  superposition  of  the  modes.  The 
POD  modes  in  this  study  are  computed  using  Sirovich's  method  of  snapshots  [19,20] 

The  main  advantage  of  the  POD  is  that  the  flow  can  be  reconstructed  using  a  selected 
number  of  orthogonal  modes  that  represent  a  progressively  smaller  portion  of  the  overall 
variance.  The  first  mode  is  equivalent  to  the  best  least-squares  approximation  to  the  wavefront 
variations,  the  second  mode  is  the  second-best  representation,  and  so  on.  POD  modes  can  be 
used  to  form  reduced-order  models  of  the  flow  by  projecting  the  governing  equations  onto  them 
[18].  They  are  also  potentially  useful  in  adaptive  optics  applications,  since  they  provide  lenslet 
array  shapes  that  would  counteract  the  largest  aero-optical  effects  [21].  In  the  present  case,  the 
POD  modes  of  the  wavefront  in  the  beam  serve  two  purposes.  First,  the  POD  modes  of  OPD 
from  the  LES/RANS  simulation  can  be  directly  compared  to  those  for  the  experimental  data, 
providing  a  statistical  comparison  of  measured  and  computed  aero-optical  aberrations.  Second, 
the  POD  modes  reveal  the  size  and  distribution  of  coherent  structure  in  the  aperture,  essentially 
identifying  turbulent  structures  in  the  flow.  However,  the  POD  does  not  produce  information 
about  the  actual  flow  dynamics,  since  each  POD  mode  contains  information  over  the  full 
temporal  spectrum. 

Figure  9  shows  the  first  10  POD  modes  of  the  mean,  tip,  and  tilt-removed  OPD  from  the 
experiment  and  simulation.  The  first  POD  mode  in  the  experimental  data  is  not  well-captured  by 
the  simulation.  However,  there  is  excellent  agreement  in  the  subsequent  modes  when  the 
simulated  modes  are  ‘shifted’  by  one  so  that  simulated  mode  1  corresponds  with  experimental 
mode  2,  etc.  The  correlation  is  flipped  somewhat  in  modes  6-8,  with  simulated  mode  6  matching 
closely  with  experiment  mode  8,  6  with  5,  and  8  with  7.  Note  that  the  modes  are  ordered  by 
eigenvalue  magnitude,  and  some  mode  swapping  can  be  expected  if  the  eigenvalues 
(representing  total  variance  contribution  of  the  mode)  are  close  enough.  In  the  simulated 
decomposition,  all  three  of  the  modes  in  question  each  account  for  approximately  5%  of  the  total 
variation  in  OPD. 

It  is  likely  that  the  discrepancy  between  the  experimental  and  simulated  optical 
aberrations  is  due  to  the  peculiarity  of  the  wall  boundary  layer  in  the  TGF.  It  is  noted  in  Smith,  et 
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(a)  Experiment  data 


(b)  LES/RANS  data 


Figure  9.  First  10  spatial  POD  modes  of  optical  path  difference  through  cavity  flow  and  wind 

tunnel  wall  boundary  layer 

al  [7]  that  the  boundary  layer  along  the  wall  in  question  is  tripped  somewhere  upstream  of  the 
test  section  and  therefore  might  display  some  wake-like  behavior  that  would  influence  the  optics. 

A  comparison  of  the  experimental  and  simulated  decomposed  OPD  modes  for  the 
isolated  boundary  layer  is  given  in  Figure  10.  Modes  1  through  3  and  6  through  10  show  very 


(c)  Experiment  data 


(d)  LES/RANS  data 


Figure  10.  First  10  spatial  POD  modes  of  optical  path  difference  through  wind  tunnel  wall 

boundary  layer 

good  agreement  with  the  experiment,  with  the  caveat  that  the  ordering  of  the  7th  and  8th  modes 
in  the  simulation  is  reversed.  The  fourth  and  fifth  modes  of  the  experiment  are  not  captured  by 
the  simulation  and  are  likely  unique  to  the  semi-tripped  boundary  layer  found  in  the  experimental 
facility.  In  particular,  the  shape  of  mode  4  of  the  experimental  boundary  layer  data  is  very  similar 
to  the  first  mode  of  the  cavity  data  and  could  reflect  the  same  turbulent  structure.  It  is  unclear 
why  this  feature  appears  to  be  more  prominent  in  the  cavity  flow  than  the  boundary  layer  flow, 
or  why  the  LES/RANS  simulation  is  better  able  to  capture  the  first  POD  mode  of  the  tunnel  wall 
boundary  layer.  It  is  possible  that  the  tripped  boundary  layer  in  the  wind  tunnel  forms  a  trapped- 
mode  resonance  when  the  cavity  apparatus  is  placed  in  the  tunnel;  this  could  potentially  amplify 
any  errant  structure  in  the  flow.  The  deviation  could  also  be  due  to  a  free-stream  disturbance  in 
the  wind  tunnel  that  propagates  into  the  interrogation  region. 

It  is  interesting  to  note  the  apparent  similarity  between  the  OPD  modes  for  the  boundary 
layer  and  cavity  flows.  In  Figure  11,  the  first  ten  OPD  modes  from  the  two  LES/RANS 
simulations  are  presented  together  for  comparison,  along  with  the  modes  of  OPD  from  the 
isolated  cavity  flow  (with  the  wall  boundary  layer  omitted  from  the  beam  path).  Taking  mode 
reordering  into  account,  nine  out  of  ten  modes  of  the  boundary  layer  and  isolated  cavity  flow 
display  a  similar  structure,  with  mode  5  being  the  exception.  For  better  clarity  of  comparison,  the 
color  map  of  some  of  the  boundary  layer  modes  has  been  reversed  -  because  the  modes  are  basis 
functions,  their  sign  is  arbitrary.  The  comparison  indicates  that  a  main  driving  factor  behind  the 
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optical  aberrations  in  both  flows  is  the  fine  eddy  structure  in  the  shear  layer  and  boundary  layer. 
The  filtering  procedure  to  remove  the  tip  and  tilt  components  from  the  OPD  snapshots  removes 
most  of  the  influence  of  boundary  layer  and  shear  layer  thickening,  as  well  as  some  shear  layer 
structure  which  is  larger  than  the  aperture  diameter.  Much  of  the  higher  order  fluctuation  that 
remains  is  mostly  a  function  of  the  free-stream  Mach  and  Reynolds  number,  which  explains  the 
agreement  in  the  higher  boundary  layer  and  cavity  modes.  On  the  other  hand,  the  large  aero- 
optical  structures  due  to  the  cavity  oscillation  dynamics  are  at  least  somewhat  reflected  in  the 
cavity  data  and  likely  cause  the  marked  deviation  of  the  first  three  POD  modes.  A  plot  of  the 
cumulative  modal  variance  contribution  (Figure  12)  shows  that  the  boundary  layer  POD  is  less 
compact,  indicating  that  more  of  the  variation  in  the  cavity  flow  data  is  due  to  the  lower,  larger 
spatial  modes. 


Combined  Boundary  Layer  and  Cavity  flow 


Cavity  flow  only 


Boundary  Layer  only 

Figure  11.  Comparison  of  first  ten  OPD  modes  of  isolated  cavity  flow  and  wind  tunnel  wall 

boundary  layer  flow 


Figure  12.  Cumulative  modal  variances  as  a  fraction  of  overall  wavefront  variance 
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4.5  Contribution  of  Various  Flow  Regions 

It  is  of  interest  to  determine  the  aero-optical  contributions  of  different  turbulent  regions 
within  the  flow,  since  this  is  information  that  cannot  easily  be  determined  from  the  path- 
integrated  OPD  data  measured  in  the  experiment.  An  estimate  of  the  OPD  contribution  along  the 

beam  can  be  obtained  from  the  root  mean  square  density,  (, p 2  ~{p)2)^2 ,  which  is  plotted  along 
the  cavity  centerline  and  along  the  beam  center  in  Figure  13.  The  highest  levels  of  density 
variation  are  in  the  shear  layer,  with  a  strong  clustering  towards  the  re-attachment  point  on  the 
aft  cavity  wall.  However,  the  density  fluctuations  also  appear  significant  both  above  and  below 
the  shear  layer.  Inside  the  cavity  itself,  the  rms  density  rises  near  the  bottom  surface,  suggesting 
that  there  could  be  some  influence  from  the  boundary  layer  formed  by  reversed  flow  circulating 
upstream  from  the  shear  layer  impingement  point. 

The  turbulent  portions  of  the  flow  can  be  divided  vertically  into  three  distinct  regions:  the 
recirculation  zone  within  the  cavity,  the  shear  layer,  and  the  region  of  small  shocks  above  the 
shear  layer.  The  vertical  extent  chosen  for  these  three  regions  within  the  optical  path  is  shown  in 
Table  3.  The  shear  layer  region  was  chosen  to  be  thicker  than  the  actual  shear  layer  to  provide  a 
buffer  that  accounts  for  shear  layer  oscillation.  The  OPD  was  computed  for  the  section  of  the 
beam  passing  through  each  of  the  three  regions  and  broken  down  by  POD  as  before.  A  diagram 
of  the  beam  split  into  the  three  regions,  along  with  the  first  nine  POD  modes  from  each  region,  is 
presented  in  Figure  14.  As  expected,  the  OPD  in  the  shear  layer  shows  evidence  of  much  finer 
turbulent  structure  which  shows  up  as  noise  in  the  POD  modes,  indicative  of  the  small-scale 
eddies  generated  in  the  layer.  The  two  most  energetic  modes  of  the  shock  and  shear  layer  regions 
are  similar  in  shape  and  determine  the  first  modes  of  the  full  line  of  sight.  The  cavity  region  also 
exhibits  these  modes,  but  in  reversed  order.  A  correlation  between  modal  shapes  is  expected 
since 


Figure  13.  Normalized  density  variance  from  LES/RANS  simulation.  Left:  Centerplane  of  cavity 

region.  Right:  Centerline  of  optical  path 
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Full  Line  of  Sight 


Cavity  Region 


Figure  14.  First  9  OPD  modes  for  different  flow  regions  within  optical  path.  Density  snapshot  along 

beam  shown  on  left. 


Table  3.  Extent  of  distinct  flow  regions  within  optical  path;  cavity  lip  located  at  y=0 


Region 

ymin  (m) 

ymax  (m) 

Cavity 

-0.037 

-0.013 

Shear  Layer 

-0.013 

0.01 

Shock 

0.01 

0.1 

Figure  15.  Cumulative  modal  variances  as  a  fraction  of  overall  wavefront  variance  for  different 

flow  regions 


Approved  for  public  release;  distribution  unlimited. 

17 


Table  4:  Computed  OPDrms,s  values  for  different  flow  regions  within  optical  path 


Region 

OPDrms, s  (nm) 

Full  line  of  sight 

14.19 

Shock  Region 

2.97 

Shear  layer 

10.34 

Cavity  region 

3.10 

the  shear  layer  is  ultimately  the  driver  of  density  variations  in  the  other  two  regions.  The  small 
shock  waves  in  the  shock  region  are  generated  by  the  mean  flow  encountering  large  eddies  on 
the  upper  shear  layer  surface,  while  part  of  the  cavity  region  turbulence  is  generated  by  air 
entrained  by  those  same  downstream-traveling  vortices.  The  structure  of  the  last  three  modes  in 
each  region  begins  to  show  large  differences;  this  is  perhaps  indicative  of  the  devolution  of  large 
eddies  in  the  shear  layer  into  smaller  structures,  which  does  not  occur  (to  this  degree)  in  the  other 
two  regions. 

The  modal  variance  distribution  (Figure  15)  indicates  that  the  shock  region  POD  is  by  far 
the  most  compact,  with  the  first  10  modes  contributing  85%  of  the  total  variance.  On  the  other 
hand,  the  shear  layer  variance  curve  is  similar  to  that  of  the  wall  boundary  layer,  achieving  90% 
only  at  the  first  60  modes.  While  the  wavefront  distortions  in  the  shear  layer  reflect  a  broad  range 
of  turbulent  spatial  scales,  the  shock  region  exhibits  regular,  largely  coherent  structure.  Despite 
the  fact  that  the  total  OPDrms  contribution  of  the  shock  region  (Table  4)  is  one  third  of  that  of 
the  shear  layer,  the  first  two  spatial  POD  modes  appear  to  be  the  primary  influence  on  the  POD 
for  the  full  line  of  sight.  This  also  means  that  any  discrepancy  between  the  modeled  and 
experimental  free-stream  conditions,  such  as  high  levels  of  free-stream  turbulence,  could 
significantly  affect  the  prediction  of  wavefront  POD  modes  if  it  disrupts  the  shock  structure.  It  is 
difficult  to  draw  hard  conclusions  only  from  the  spatial  modes.  However,  it  is  clear  that  the 
contribution  from  the  shock  and  cavity  zones  is  important  to  the  overall  aero-optics  of  this 
problem,  and  that  even  after  removing  the  jitter  caused  by  bulk  shear  layer  oscillation,  modeling 
this  problem  as  a  simple  mixing  layer  is  insufficient  from  the  aero-optical  point  of  view. 


4.6  Dynamic  Mode  Decomposition  of  OPD 

POD  analysis  can  be  useful  as  a  means  of  data  comparison  and  for  identifying  spatial 
modes,  but  does  not  necessarily  give  insight  into  the  dynamics  of  a  system.  Several  recent  aero- 
optics  studies  [20]  have  used  Dynamic  Mode  Decomposition  (DMD)  to  analyze  the  spectral 
behavior  of  flows.  Unlike  POD,  DMD  provides  a  modal  breakdown  of  the  data  where  each  basis 
function  is  tied  to  a  single  frequency.  This  is  of  obvious  use  in  the  analysis  of  systems  driven  by 
discrete  forcing  frequencies,  such  as  the  Rossiter  modes  in  a  cavity  flow.  In  theory,  DMD  would 
allow  the  isolation  of  flow  dynamics  that  occur  at  the  key  frequencies. 

The  standard  DMD  algorithm  as  given  in  Hemati,  et  al.  [22]  is  as  follows:  as  with  POD,  we 
compile  a  matrix  of  snapshot  vectors  X  =  [x1,x2,...,xJV_1].  A  second,  offset  snapshot  matrix  is 
formed  such  that  Y  =  [x2,x3,...,xAr].  The  first  snapshot  matrix  undergoes  singular  value 
decomposition, 


X=QxXV' 


(20) 
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After  which,  the  DMD  operator  matrix  is  constructed: 


K  =  Q'xYX+Qx  =Q'xYY2_1  (21) 

where  X+ is  the  Moore-Penrose  pseudoinverse  of  X.  The  DMD  modes  and  frequencies  are  then 
given  by  the  eigenvectors  and  eigenvalues  of  K ,  such  that  if  Vj  is  the  jth  eigenvector  with 
eigenvalue  Xj Then  Qxv;is  the  jth  DMD  mode  with  frequency  /•  =  imag(]og(Aj )) /  2nAt  with  At 

being  the  time  interval  between  successive  snapshots. 

In  this  work,  an  alternative  ’streaming1  DMD  algorithm  developed  by  Hemati,  et.  al  [22] 
is  used.  Instead  of  constructing  the  snapshot  matrices  in  a  single  pass,  the  algorithm  appends  a 
single  pair  of  snapshots  to  the  system  at  a  time,  updating  the  basis  sets  Qx  and  the  corresponding 
Qr  as  snapshots  are  added.  Limiting  the  ranks  of  the  basis  sets  enables  POD  compression  of  the 
data  set,  which  reduces  noise  in  the  snapshot  data  and  allows  the  algorithm  to  pick  out  a  small 
number  (equal  to  the  specified  rank)  of  ’key’  modes  important  to  the  system  dynamics,  rather 
than  a  full  spectrum  of  modes.  This  is  advantageous  in  the  current  problem  because  the  limited 
time  frame  of  the  data  collection  and  high  levels  of  noise  makes  dominant  modes  very  difficult  to 
identify  solely  based  on  their  magnitude.  However,  it  does  introduce  uncertainty  since  the 
calculated  modes  are  sensitive  to  the  choice  of  maximum  rank.  Presently,  the  maximum  rank  of 
the  DMD  operator  was  set  to  50,  which  was  sufficient  to  isolate  all  of  the  larger  peaks  in  the 
OPD  spectra  along  the  optical  beam  path. 

The  DMD  was  applied  to  snapshots  of  OPD  spanning  the  entire  cavity,  which  were 
sampled  at  24  kHz.  The  total  sample  size  was  215  snapshots.  The  choice  of  the  larger  aperture 
size  serves  to  better  give  an  idea  of  the  overall  flow  dynamics  which  govern  the  OPD  over  the 
cavity,  rather  than  a  small  sample  at  the  cavity  center.  Furthermore,  it  is  not  guaranteed  that  the 
3cm  aperture  of  the  beam  is  large  enough  to  capture  the  largest  aero-optical  structures  in  the 
flow.  As  before,  the  flow  was  split  into  three  distinct  vertical  regions  spanning  the  recirculation 
zone  in  the  cavity,  the  shear  layer,  and  the  shock  region  above  the  cavity.  Figure  16  compares  the 
spectra  of  the  computed  DMD  modes  with  an  average  OPD  power  spectral  density  over  the 
optical  beam  path.  It  should  be  noted  that  the  OPD  spectra  over  the  full  cavity  do  not  necessarily 
match  those  over  the  smaller  aperture;  the  latter  were  used  mostly  because  they  were  sampled  at 
a  higher  frequency  of  106  kHz  and  thus  confirmed  that  the  detected  peaks  above  10  kHz  were 
not  solely  due  to  aliasing. 

The  shock  region  spectra  (Figure  16a)  show  distinct  peaks  at  the  third,  fourth,  and  fifth 
Rossi  ter  modes  (with  the  third  mode  being  dominant),  as  well  as  a  peak  around  11  kHz  before  a 
rapid  drop-off  in  intensity  occurs  at  higher  frequencies.  The  dominant  modes  computed  by  the 
streaming  DMD  algorithms  cluster  near  the  Rossiter  modes  and  the  1 1  kHz  peak.  The  shear  layer 
spectra  (Figure  16b)  displays  a  much  more  uniform  distribution  of  energy,  with  a  higher  relative 
PSD  at  frequencies  into  the  tens  of  kilohertz.  Relative  peak  intensities  are  seen  at  the  first  five 
Rossiter  modes.  The  DMD  produced  peaks  close  to  the  Rossiter  modes  as  well  as  a  largely 
uniform  clustering  between  5.5  and  11  kHz.  In  contrast,  the  recirculation  region  spectra  (Fig. 
16c)  are  dominated  by  the  lower  frequencies,  with  dominant  peaks  at  the  first  two  Rossiter 
modes  as  well  as  a  frequency  in  between. 

The  real  components  of  the  main  DMD  modes  for  each  of  the  flow  regions  are  presented 
in  Figures  17-19.  The  black  circle  on  each  mode  shows  the  location  of  the  optical  beam  for 
reference.  Overall,  the  shock  region  exhibits  the  smoothest  and  most  coherent  optical  structures 
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in  the  form  of  regular  waves  spanning  the  cavity  which  increase  in  magnitude  as  they  travel 
downstream  and  decrease  in  wavelength  at  higher  frequencies.  These  mode  shapes  roughly 
correspond  with  the  shear  layer  modes  (Figure  18)  at  similar  frequencies,  although  the  latter  are  a 
lot  noisier.  The  structural  similarity  is  consistent  with  the  formation  process  of  the  shocklets 
above  the  cavity,  which  occur  when  the  shear  layer  deflects  into  the  free-stream.  Thus  the  aero- 
optical  structures  in  the  shock  region  assume  the  shape  of  the  larger-scale  oscillations  of  the 
shear  layer  and  propagate  them  into  the  free-stream  above  the  cavity.  In  essence,  the  shock 
region  structures  are  an  inviscid  footprint  of  the  ones  in  the  shear  layer.  This  relationship  could 
factor  into  the  design  of  an  effective  adaptive  optics  system  for  supersonic  cavity  flows  -  namely, 
a  configuration  which  reduces  the  aero-optical  aberrations  due  to  shear  layer  oscillation  at  the 
Rossiter  frequencies  would  at  the  same  time  potentially  counteract  the  aberrations  caused  by  the 
shock  region. 

The  dominant  modes  in  the  recirculation  region  (Figure  19),  on  the  other  hand,  show  less 
uniform  behavior.  The  lowest-frequency  mode  appears  almost  like  a  standing  wave  spanning  the 
entire  length  of  the  cavity.  Near  the  second,  third,  and  fourth,  Rossiter  modes  (f2,  f3,  and  ft)  there 
is  an  apparent  split  at  the  cavity  midpoint.  In  the  downstream  half  of  the  cavity,  the  modes  show 
large  structures  with  high  intensity  that  do  not  span  the  cavity  width  and  are  due  either  to  the 
impingement  of  the  shear  layer  on  the  rear  cavity  wall  or  vortices  caused  by  large  structures  in 
the  shear  layer  entraining  air  from  the  recirculation  region.  In  the  upstream  half  of  the  cavity,  the 
mode  shapes  resolve  into  relatively  faint  bar  structures  indicative  of  upstream-traveling  pressure 
waves.  At  the  highest  detected  frequency  (f5  =  9.390  kHz),  the  DMD  mode  consists  of  isolated 
structures  that  shrink  downstream.  The  physical  process  behind  these  is  not  immediately  obvious 
but  could  be  due  to  the  influence  of  convecting  structures  in  the  shear  layer. 

5.0  CONCLUSIONS 

In  this  study,  a  hybrid  LES/RANS  methodology  was  applied  to  a  high  Reynolds  number 
supersonic  open  cavity  flow  with  an  optical  aperture,  with  the  driving  motivation  being  to 
confirm  that  the  low-dissipation  LES/RANS  approach  can  provide  accurate  simulation  of  the 
small-scale,  transient,  turbulent  effects  that  have  a  large  impact  on  the  performance  of  optical 
systems.  The  simulation  showed  good  agreement  with  Schlieren  photography  from  the 
experiment  and  the  shear  layer  oscillation  was  well-predicted  by  Rossiter's  empirical  formula. 
The  root  mean  square  OPD  was  over-predicted  by  approximately  10%  in  the  simulation.  The 
proper  orthogonal  modes  of  the  wavefront  distortion  compared  very  well  with  the  experimental 
data  with  the  exception  of  a  prominent  spatial  structure  seen  in  the  experiment,  which  may  be 
due  to  wind  tunnel  noise.  Despite  prior  studies7  suggesting  anomalous  behavior  in 
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(b)  Shear  layer  region 


Figure  16.  Comparison  of  DMD  mode  frequencies  and  aperture-averaged  OPD  spectra 
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fi=  0.883  kHz  f2=  1.988  kHz  f3  =  3.214  kHz 


f4=  10.947  kHz  f5  =  5.713  kHz  f6=  10.652  kHz 

Figure  17.  Full-cavity  shock  region  DMD  modes 


fi  =  0.717  kHz 


U  =  3.845  kHz 


f5  =  5.236  kHz 


f3  =  3.320  kHz 


f6=  10.280  kHz 


Figure  18.  Full-cavity  shear  layer  region  DMD  modes 
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f  i  =  0.750  kHz  f2  =  1 .479  kHz  f 3  =  1 .949  kHz 


U  =  3.324  kHz  f5  =  4.906  kHz  f6  =  9.390  kHz 


Figure  19.  Full-cavity  recirculation  region  DMD  modes 

the  boundary  layer  over  the  wind  tunnel  window,  computed  POD  modes  and  rms  OPD  of  the 
simulated  boundary  layer  were  in  close  agreement.  Overall,  the  results  suggest  that  an 
LES/RANS  approach  is  sufficiently  robust  to  accurately  model  aero-optical  aberrations,  even  for 
non-idealized  flows. 

A  further  decomposition  of  the  cavity  flow  into  distinct  regions  showed  that  while  the 
shear  layer  spanning  the  cavity  was  responsible  for  the  bulk  of  the  aero-optical  distortions,  the 
region  of  shock  waves  caused  by  shear  layer  entrainment  into  the  free-stream  also  produced  a 
significant  contribution,  and  was  responsible  for  some  of  the  coherent  structure  seen  in  the  POD 
of  the  wavefront.  A  dynamic  mode  decomposition  of  a  much  larger  wavefront  spanning  the 
entire  cavity  was  performed  and  showed  that  the  main  driver  of  aero-optical  aberrations  is  the 
bulk  flapping  motion  of  the  shear  layer  occurring  at  the  Rossiter  frequencies.  The  oscillation 
causes  a  similarly- shaped  wavefront  distortion  within  the  shock  layer  region  over  the  cavity.  At 
lower  frequencies,  the  structures  larger  than  the  aperture  of  the  Shack-Hartmann  sensor  used  to 
measure  OPD  in  the  experiment,  so  a  significant  aperture-size  effect  is  to  be  expected. 

There  are  several  uncertainties  left  to  explore  regarding  the  use  of  LES/RANS  models  to 
simulate  aero-optics  problems.  Foremost,  while  it  is  assumed  that  the  significant  aero-optical 
aberrations  occur  in  the  free-stream,  the  actual  effect  of  the  near-wall  RANS  model  choice  and 
blending  function  has  not  been  studied.  A  related  problem  is  with  the  modeling  of  flow  near  the 
wall  with  the  immersed  boundary  method.  Although  in  the  present  case  the  assumption  that  the 
flow  developing  upstream  of  the  cavity  is  fully  turbulent  appeared  to  produce  adequate  results, 
this  is  far  from  a  safe  assumption.  For  example,  if  the  point  of  interest  was  aero-optical 
distortions  over  the  flat  plate  upstream  of  the  cavity,  the  precise  IB  wall  model  would  certainly 
be  consequential.  Nevertheless,  if  the  near-wall  modeling  is  handled  in  a  careful  way,  the  IB 
methodology  coupled  with  LES/RANS  is  potentially  attractive  for  simulating  flows  over  more 
complicated  geometries,  such  as  optical  turrets. 
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